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ABSTRACT 

A way to probe alternative theories of gravitation is to study if they could 
account for the structures of the universe. We then modified the well-known 
Gadget-2 code to probe alternative theories of gravitation through galactic 
dynamics. As an application, we simulate the evolution of spiral galaxies to 
probe alternative theories of gravitation whose weak field limits have a Yukawa- 
like gravitational potential. These simulations show that galactic dynamics 
can be used to constrain the parameters associated with alternative theories of 
gravitation. It is worth stressing that the recipe given in the present study can be 
applied to any other alternative theory of gravitation in which the superposition 
principle is valid. 

Subject headings: galaxies: spiral - gravitation - methods: numerical 



Introduction 



One of the most important challenges of modern cosmology concerns the dark energy 
problem. The nature or origin of dark energy cannot be associated with particles, and 
i ts interpretation b ased on quantum field theories does not provide a satisfactory answer 
( ICarroU et al.lll992l ). The scalar field hypotheses, as given by quintessence models, do not 
either solve the problem satisfactorily, but instead, give rise to other questions without 
explanations. Moreover, observations claim that ~ 70% of the total energy composition of 
the universe is r nade of this puzzling (c osmological) ingredient that accelerates the expansion 
of the universe (jPerlmutter et al.lll999l ). like an antigravity term in the Einstein's equations. 
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Another puzzling ingredient of the Cosmos is the dark matter ( iBahcall et aLl |2004| ). 
Cosmology shows us, based on observations of galactic and cluster dynamics, gravitational 
lenses, etc., that the dark matter composition of the universe is ~ 25%. In principle, the 
dark matter nature can be explained by particle physics; one expects that the large hadron 
collider helps us to give a good answer to this issue. 

Recently, there appear in the literature some new alternative theories to the Einstein's 
General Relativity aiming to explain the structures of the universe without dark contents 
(matter and energy). They are based on different hypotheses (e.g., massive gravitons. 



I (see, e.g.. 


VIoffat & Sokolov 


1996 




Piazza & Marinoni 


Sienore 


2005 


de Araujo & Miranda 


2007 


)• 



In some theories, for example, the gravitational potential is Yukawa-like (hereafter 
Yukawian gravitational potential, YGP) in the weak field limit. The YGP reads 



Gm 



~r/A 



where m is the point-mass source of the field, r is the distance to the point mass, and A is a 
characteristic length. In some theories, one postulates the existence of a massive boson called 
graviton, of mass mg, whose Compton wavelength can be interpreted as the A parameter. 

It is worth noting that many authors consider that to probe a given alternative 
potential it is enough to reproduce the rotation curves of spiral galaxies assuming centrifugal 
equilibrium. In fact, the best way to probe a potential is to use a galactic dynamical 
approach. Moreover, it is necessary to verify, even using simulations (living systems), not 
only the rotation curves but also the surface density profiles. From the observations one 
can infer that this profile is exponential, therefore a given alternative potential must be 
consistent with it too. 

Note also that alternative gravitational potentials are basically of two types: those 
in which the superposition principle applies and those in which it does not. For the first 
type, one can use iV-body simulations, in which the use of the superposition principle is 
inherent. For the second type one must use, for example, smoothed particle hydrodynamics 
simulations, where it is not necessary to be concerned with the superposition principle. 
Modified Newtonian dynamics, for instance, is of the second type, therefore one cannot use 
A^-body simulations to model mondian galaxies. 

In our previous papers ( Brandao &: de Araujo 2010al lbl). we modified and tested an 
A^-body code replacing the Newtonian potential by the YGP. Moreover, we made some 
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numerical simulations to study elliptical galaxies. It is worth mentioning that although we 
adopted this particular potential, the recipe given in these papers can be applied to any 
alternative gravitational potential in which the superposition principle is valid. 

To complete our studies of galactic systems with alternative gravitational potentials, 
we consider in the present paper the modeling of late-type systems. Therefore, our main 
aim here is to give a recipe to perform such a study, and we adopted again the YGP as an 
example. 

Concerning in particular the YGP, in almost all previous works concerning it, the 
investigations had be e n ma de using analytical or numerical approaches. For example. 



de Araujo fc Mirandal (120071 ) have probed how the YGP can change the rotation curves 
of spiral galaxies, using "static" (centrifugal equilibrium) and analytic models of exponential 
disks. 

In the present work, as an alternative to the approach used by de Araujo & Miranda, we 
study numerical A^-body simulations of disk galaxies and investigate how the YGP changes 
the canonical morphology of a simulated disk. 

We show in the end how "living" systems behave under the YGP. Verifying therefore 
whether the YGP can or cannot keep the disk galaxies in dynamical equilibrium. 

The paper is organized as follows: in Section 2, we present the galactic model used 
and the numerical code adopted to perform the simulations; in Section 3, we present the 
simulations and discuss the results; and, finally, in Section 4 we present the main conclusions. 



METHOD AND SCENARIO: A SPIRAL GALAXY MODEL UNDER 

YGP 



To investigate YGP at galactic scale, or any other alternative potential in which the 
superposition principle can be applied, one has to choose a typical model (dynamical 
equilibrium) for the galaxies, and write an efficient A^-body code, based on the tree method 
jBarnes fc Hutlll986h : or example, to follow the evolution of the galaxies. 



In the present paper we will consider disk galaxies. We know that such galaxies, 
their observational, structural and dynamical properties are well understood by the galactic 



dynamics approach (see, e.g.. iBinney fc Merrifieldlll998t iBinney fc Tremainell2008l ). These 



systems are well modeled by numerical simulation tools and their secular evolution, under 
the mutual forces arnqng their particles, can be followed (se e, ag^ 



Springel fc Whitelll999l ISpringel et al.ll2005l iRomero-Gomez et al.ll2006 



Hernguist 



Athanassoula 



1993 



2006 
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and references therein). 

In the preseii t pape r, our model of disk galaxies is based on the model s described by 



Springel fc White! (jl999[ ) and on some other implementations introduced by ISpringel et al. 
(l2005l ). Our galaxy model is called, throughout this paper, Springel-Di Matteo-Hernquist 
disks (hereafter SdMH disks), due to their numerical prescriptions, developed specially to 
model late-type systems. However, we here simplify the SdMH disks, since we do not include 
bulges or gas particles, due to the fact that our investigation is aimed only to study the overall 
dynamical and gravitational properties of disk galaxies under alternative potentials. 

It is worth noting again that, although we consider here the particular case of the YGP, 
the prescription considered throughout the present paper holds for any other alternative 
potential one would like to probe, since they respect the superposition principle. 



2.1. Modeling Galcixies with an A^-dody Code 



To simulate a spiral galaxy under YGP, we have chosen the Gadget-2 code (ISpringel 



20051 ) and changed its structure to include the YGP, as we show in detail in our previous 
work (IBrandao fc de Araujdl2010al ). 



Gadget-2 is based on the tree code method (IBarnes &: Hutlll986l ). So, its computational 
effort is 0(A^log(A^)) , instead of OiN"^) operation s required by direct sum algorithms. With 
this kind of code (see lBrandao &: de Arauidl2010a] . for details), we can integrate all equations 
of motion of a set of collisionless particles and foll ow their evolution, as we did for example 
for elliptical galaxies (IBrandao fc de Araujdl2010bl ). 



We recall here our previous arguments to justify our methodology. We constructed 
galaxies initially with a Newtonian potential and then we submit them to the YGP. It 
is important to emphasize that particles represent physical observable quantities, such as 
positions, velocities and masses distributed over a given volume. So, when we realize an 
initial galaxy snapshot with the Newtonian potential, we mimic the following observational 
characteristics: radial luminosity profile, radial density profile, and velocity dispersions. It 
is important to bear in mind that Newtoni an galaxies are consistent with obs ervations, if 
dark matter is taken into account (see, e.g., lOh et. alll201ll : iGuedes et. all 1201 if ). 



With this philosophy, independently on the physics used to build up galaxies, the 
simulated particles must reproduce the observed characteristics of real objects. Our aim 
is then to check, at the end of our simulations for a given alternative gravitational potential, 
if these characteristics are really consistent with observed objects. 
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One could argue that the best way to model Yukawian galaxies would be to consider 
their formations starting with collapsing halos. But since, as mentioned above, we start 
with a Newtonian model, this procedure could in principle be considered a limitation of the 
present approach. In the near future, however, we intend to investigate alternative potentials 
starting with collapsing halos. 



2.2. Assembling Galactic Halos and Disks: The particle's Positions 

Our late-type system is composed by a dark matter halo modeled by a Hernquist sphere 
and an exponential- Spitzer disk. The halo density distribution law reads 

Mdm a 

Pdm(r) = — , (2) 

ZTT r[r + a)'^ 

where pdm{^) is the radial density profile of the dark matter (dm) halo, Mdm is the halo 
mass, r is the radial distance from the center of the whole mass distribution, and a is 
a characteristic length of the halo's core. This parameter is related t o the concentration 



param eter c of the Navarro, Frenk and White (hereafter NFW) sphere (INavarro et al.lll996 



19971 ) of concentration index c = r2oo/fs, where r2oo is the virial radius, where the mean 
overdensity, as compared to the critical density, is 200. The Vs parameter is a characteristic 
scale length of the NFW sphere. 



as 



Recall that the scale lengths of the Hernquist sphere and the NFW sphere are related 

a = r, V2[ln(l + c)-c/(l + c)]. (3) 



Also, recall that the mass of the NFW sphere diverges as r — )■ oo, while the mass of the 
Hernquist spheres converges as r — )■ oo to the value Mdm- The two profiles agree very well 
within the virial radius r2oo- But NFW spheres mu st be truncated, while Hernquist spheres 



do not. In this way, we follow ISpringel et al.l (120051 ) and use the Hernquist sphere to model 
the dark matter halo. 

In the present work, we consider a disk made of baryonic particles. Our model does not 
consider gas, star formation, wave shocks, supernovae or supermassive black holes. We have 
also considered that the baryonic mass is a fraction of the total mass Mjisk = ""^dAftot, 
where rrid is a dimensionless parameter, Mdisk is the disk mass, and the total mass is 
Mtot = Mdisk + Md^. 

The disk has a residual spin from its primordial halo, from which it was formed. The 
spin parameter reads 
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A. 



(4) 



where J is the total angular momentum of the primordial halo, E is the total energy, G is 
the universal gravitational constant, and M = Mdm + ^disk is the primordial mass. The spin 
parameter is used to compute the parameter fc that will be defined below. 

The disk density distribution is given by an exponential law and the Spitzer isothermal 
sheet, namely: 

where Zq is the disk t hickness, h is the r adial scale length, and is the total disk 
mass. Differently f rom ISpringel et al.l (120051 ). that left Zq as a free parameter, we followed 
Springel fc White! ( 119991 ) and set zq ~ 0.2Rd due to the fact that many late-type systems 
have this typical scale length. All particle's positions are distributed by the Monte Carlo 
method. 



2.3. Assembling Galactic Halos and Hisks. The Particle's Velocities 



The main ingredient of this subsection has to do with how to distribute velocities of the 
halo and disk particles. While the entire prescription is found in the literature, we recall 
here only the most important ingredients to model our galaxies. 

The velocity structure of our model depends on the calculation of the potential generated 
by the matter. Hernquist spheres have a potential given by 



CM, 



dm 



dm 



(6) 



To calcula t e the potential from the disk contrar y to ISpringel fc White (|l999[) and 
Springel et al. ( 2005 ). we follow Equation (2.170) from Binney fc Tremainel ( 2008 ). who 
consider an exponential thick disk of a completely flattened homoeoid. 

The calculation of the potentials is a key issue, since it is used to compute the velocity 
dispersions. For this axisymmetric system, we as sume that the velocity distribution function 



is f{E, Lz) (see, e.g., iMagorrian fc Binneylll994l ). where E is the total energy and is the 



^-component of the angular momentum. 



It follows that the flrst velocity moments are given by (in cylindrical coordinates) 
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VR = Vz= VrV, = V,V^ = VrV^ = 0, (7) 



V 



2 



4 = ^z, 



+ (10) 



where the circular velocity is Vc = -Rf^; with the bars denoting the mean over the quantities 
in consideration. In the above equations, p is the density of the quantity for which we 
compute the velocity variances, and the potential is due to the whole matter distribution. 

From the distribution function /(-R, L^), we conclude that in the azim uthal direction the 



mean streaming velocity in not necessary null. We, therefore, follow ISpringel fc White 



(119991 ) ■ considering = fgVc, where generally /s(^ 1) is a factor that depends on Ag. This 
means that the streaming velocity of the dark matter halo is a fraction of the local circular 
velocity. Once specified all the above values, the velocity dispersions for the dark matter 



halo are given by cTj = \/ (vf), with i = z,R, and al — _ ,2 



V2 — 



To the disk, the calculations are similar to the presented above, although the calculation 
of the component is very different. First, the mean streaming is estimated by the epicyclic 
approximation, and the equations used to calculate the velocity dispersions are given by 



2 

f^^ = — , (11 



where 



2 



4 9$ /3 9$ 



^ RdR \ RdR 



Once these quantities are obtained, we use Equation [TO] to evaluate the streaming 
velocity: 
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a 



R 



7]^ 



1/2 



(13) 



In a few words, all this prescription consists in calculating ak, with k = R,z,(f) for the 
disk and for the halo at any {R, z) points. To this aim, we employ the following computational 
techniques. 



We build a logarithmic mesh, where the density and potentials are calculated at the 
respective (i?, z) points. 

We use subroutines based on the spline techniques to make the above integrals. 
We interpolate the dispersions at the particle's points. 

The particle velocities are s et by random numbers from the Schwarzschild's distribution 
(IBinney fc Tremaindl2008l ). which is given by 



f{v)d^v 



NdH 



exp 



(27r)3/V^a^a 
where is the number of particles per unit volume. 



+ 



+ 



2a| 2al ' 2a? 



(14) 



Following the above prescription, we have written a code to build up a dark matter halo 
and a baryonic disk particle. We consider both kind of particles gravitationally coupled. We 
set iVhaio = 30,000 particles for the dark matter halo and iVdisk = 30,000 for the disk in our 
models. We will see later on that a higher resolution was also used in some simulations. 



^The following set of default parameters are chosen to realize SdMH disks ( ISpringel et al. 

2005h : total mass Mt = v^Qf^/ilOG Ho) = 0.98 x lO^^^^^o, where waoo = 160A;m s"^ is the virial 
velocity, G is the gravitational Universal constant, Hq = lOOkm Mpc~^ is the Hubble 
constant; the total mass of the disk = mdMt, where rrid = 0.041 is a dimensionless 
fraction of the total mass, the disk scale length h = 2.74 kpc, the disk vertical scale-height 
Zq ~ 0.2/1, and the spin parameter A = 0.033. 

In Figure [H we display the rotation curves of the modeled galaxy , using the parameter s 
described above. We note that our results are very similar to that by ISpringel et al.l ( l2005l ). 
The differences reside on the numerical procedures and techniques to compute the disk 
potentials, particle noise due to the halo, disk truncation, etc. 
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3. SIMULATIONS AND DISCUSSION OF THE RESULTS 

We now present and discuss the results of our simulations. As usual in treecodes, we 
have chosen the tolerance parameter 9 = 0.8, which gives a better performance for the 
calculations. Other typical parameters are the halo smoothing scale length Ih = 0.15 and 
the disk smoothing scale length 1^ = 0.10. It is worth noting that we have changed the scale 
length parameters in the following range: Ih = 0.15 ± 0.1 and Id = 0.10 ± 0.05; and our 
results do not change significantly. 

The following simulations are performed: Newtonian potential with the default Gadget- 
2 code and YGP with A = 1, 10, 100, and 1000 kpc. 

The justification for choosing the above values for A is the following. Since spiral galaxies 
have characteristic dimensions of tens of kpc we choose the Yukawa scale lengths to be much 
lower (i.e., 1 kpc), much greater (i.e., 1000 kpc), and of the order (i.e., 10 and 100 kpc) of the 
spiral galaxy sizes, in order to see how the galaxy structure is affected. It is expected that 
for A much larger than the characteristic dimensions of galaxies, the YGP model is similar 
to the Newtonian one. On the other hand, for A smaller than the characteristic dimensions 
of galaxies, the YGP and the Newtonian models yield very different results. 

3.1. The Newtonian Simulation 

We use the default Gadget-2 code and make a typical Newtonian simulation with a 
total simulated time of t = 1 Gyr. This model is used like a "control group" , with which we 
compare the other runs. 

In Figure [21 we show the principal features of the Newtonian run. In this figure we 
display at the top left panel the phase space points r x f of the initial snapshot, where r is 
the distance from the center of the matter distribution in kpc and v is the modulus of the 
velocity, in kms~^; at top right, we show the relative energy conservation AE/Eq, where 
E{t) is the total energy of the simulated system (disk plus halo) at time t, AE = E{t)—E{0), 
and E{0) = Eq. From this frame, we conclude that the energy conservation violation is less 
than 1%, showing that our simulation is reliable. At the bottom left, we present the phase 
space points at final time, t = 1 Gyr; Finally, at bottom right, we present the rotation curve 
Rx Vr, where R is the radial cylindrical coordinate and Vr is the rotation velocity. Note the 
similarity between the initial and the 1 Gyr rotation curves; one concludes that early-type 
galaxy structure can be accounted for the Newtonian model quite well, as is well known. 

In Figure [3l we show the snapshots at t = 0, 0.33, 0.66 and 1 Gyr. We note that the disk 
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evolves to a late-type system with a central bar, for which spiral arms are present in each 
end. Bars usually appear in A^-body simulations, and their developments in these simulations 
seem to be connected with the various parameters listed in Section I2.3[ It is worth noting 
that there is a branch in the Galactic Dynamics which studies the gl obal stability of the 
differ entially rotating disks and the formation of bars (see chapter 6.3 of iBinney &: Tremaine 

mi- 



lt could be argued that some characteristics of our numerical procedures could contribute 
to the bar format ion, but this feature is far from being a bad result. A recent study 
(jVerley et al.l 120071 ). based on observational data of isolated galaxies, interprets bars and 
arms as a natural consequence of secular evolution of late-type systems, and this is also 
corroborated by numerical simulations. These authors conclude that isolated galaxies do 
not seem to be preferentially barred or unbarred. Therefore, in this work, bars will be 
considered as a natural feature of the simulated late-type systems. 

In Figure m we note the develo pment of spiral arms in th e first 0.33 Gyr of simulated 
time, due to swing amplification ( iBinney fc Tremaind 120081 ). This result, as a lready 
mentioned, is expected in A^-body simulations of disk galaxies (jSpringel et al.ll2005l ). We 
conclude, from all these figures, that our Newtonian simulation yields a typical morphology 
found in the Cosmos. We emphasize that we built up other models which are stable against 
bars (mostly the big disk ones), but we chose to consider the Milky- Way like model in the 
present work. 



In the next subsections we simulate YGP late- type galaxies for different values of A. 



3.2. The YGP Simulations: A = 100 and A = 1000 kpc 

We run the galaxy model with A = 100 and A = 1000 kpc, whose results are remarkably 
similar. Therefore, we discuss here only the A = 100 kpc case. 

From Figure [5l we note its resemblance to the Newtonian simulation. In this figure, only 
the final rotation curve is slightly different from the Newtonian runs, but it is compatible 
with observational curves. In the bottom left panel of this figure, one can see a peak in the 
rotation curve around kpc < R < 2 kpc, which is due to a bar that rotates like a rigid 
body. 

The Figures [6] and [7] display the particle's positions of the disk in the xy-plane, at 
different simulated times. In particular. Figure [6] displays four snapshots for t = 0, 0.33, 0.66 
and 1 Gyr. We note that the system evolves to a disk with a central bar and some spiral 
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arms, although these arms are not so remarkable, in comparison to the Newtonian case. 

The Figure [7] shows the first 320 Myr of the disk's evolution. We can see at t = 0.08 and 
t = 0.16 Gyr some remarkable spiral arms, which are expected features, due to the swing 
amplification. We can see that the number of spiral arms changes as the simulated time 
increases, and the system evolves to a central bar earlier than in the Newtonian case. 

Although the A = 100 kpc simulation describes quite well the morphological aspects of 
spiral galaxies, it is important to know if some characteristics discussed above, such as the 
bar for example, depend on the resolution adopted (~ 10^ particles). We then construct a 
model with a higher resolution using the same physical parameters, but now with 300,000 
particles for the disk and 600,000 for the halo, respectiveljl^]. The smoothing scale lengths 
were recalculated and now read Ih = 0.01 and Id = 0.004, respectively (see Brandao & de 
Araujo 2010a, and references therein). 

In Figures [8] and m we present the results for the high- resolution simulations for A = 100 
kpc. Note that the relative energy conservation at the end of the simulation is ~ 1%, showing 
that our simulations are reliable. Comparing these figures with those for lower resolution 
one notes some similarities, although the higher resolution figures obviously stand out more 
clearly the features of the modeled galaxy. The bar still appears, showing that it is not an 
effect related to the low resolution. The rotation curves for both resolutions are also similar, 
with both resembling a Newtonian rotation curve. 

In conclusion, we have seen in this subsection some interesting results that show that 
the YGP can reliably model a spiral galaxy if A ^ 10 kpc. In this case, the YGP presents 
results that are similar to the Newtonian ones. 



3.3. The YGP Simulation: A = 10 kpc 

In Figure [TD], we display the same as in Figure E] but for A = 10 kpc. Note that the 
initial and final rotation curves are quite different. This shows that it is not possible to 
produce a A = 10 kpc Yukawian spiral galaxy consistent with the observations. 

Note that the energy violation, also shown in Figure dUl is better than in the case for 
A = 100 kpc; this implies, therefore, that our simulations are quite reliable. 

In Figure [TTl we show how would be a A = 10 kpc Yukawian spiral galaxy. In this 



^AU the higher resolution simulations were performed in the "HPC Bull Cluster" belonging to the State 
University of Santa Cruz, which was sponsored by FAPESB. 
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sequence of snapshots, we can see the growth of the disk as well as the evolution of the 
central part of the galaxy. It is interesting to note that the core becomes smaller than 
initially. The YGP for A = 10 kpc makes the central parts of the disk shrink. The first 320 
Myr of the simulation is shown in Figure [T21 We see again the swing amplification and the 
development of spiral arms. 

Similarly to what we did for the simulations described in the previous subsection, we also 
simulate a high-resolution model for A = 10 kpc. Figures [12] and [H] show the results of the 
high-resolution simulations. As in the case for A = 100 kpc, the relative energy conservation 
at the end of the simulation is ~ 1% for A = 10 kpc, showing that our simulations are 
reliable. As before, comparing these figures with those for lower resolution one notes that 
the high resolutions obviously stand out more clearly the features of the modeled galaxy, but 
the main conclusions presented above are the same. In particular, the rotation curves for 
both resolutions are similar, but they are quite different from a Newtonian rotation curve. 

Summing up a YGP with A = 10 kpc does not model a spiral galaxy appropriately. 

3.4. The YGP Simulation: A = 1 kpc 

When a spiral galaxy is submitted to the YGP, atypical morphologies appear, in 
particular for small values of A. Figures [TSlfTTI display the results for the Yukawian disk 
galaxy simulation for A = 1 kpc. 

In Figure [T51 the top right panel shows us that energy violation is very small, proving 
the reliability of this simulation. The top left panel shows that the initial positions and 
velocities of the particles in the phase space resemble those of the Newtonian simulations 
because we are using the same initial snapshot. But, when the YGP and its corresponding 
acceleration are considered, their exponential factor plays a decisive role: for A = 1 kpc, 
the forces between distant particles 1 kpc ) are almost "turned off". These particles 
become almost free and, as a result, their initial velocities are almost unchanged. As time 
goes on, the average distance between particles increases, these particles leave the galaxy 
and the system becomes more and more diffuse. Faster particles reach distant regions first, 
and the slower ones spend more time to move through the galaxy. The phase space of the 
final snapshot is displayed at the bottom left panel and shows us that the particles escape 
from the galaxy and the initial information is lost. The bottom right panel shows us that the 
initial rotation curve is lost, namely, the initial and final rotation curves are quite different. 
From these considerations, we conclude that a putative A = 1 kpc Yukawian spiral galaxy is 
ruled out, since it does not resemble the spiral galaxies observed in the universe. 
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In Figure [16] we show a face-on view of the simulated galaxy, which is gradually destroyed 
and dispersed in the intergalactic environment. See also Figure [T71 which shows what is 
happening with the central part of the disk. As time goes on, the central part becomes 
more and more empty of particles. This result leads one to conclude again that a A = 1 
kpc Yukawian spiral galaxy is ruled out. Although not shown, a high-resolution simulation 
corroborates these conclusions. 



It is worth recalling that Ide Araujo fc Miranda! ( 120071 ) showed that depending on the 
ratio between A and the scale length of the disk, the disk is destroyed. Such a result is in 
full agreement with our simulations for A = 1 kpc. 



3.5. The Surface Density Profile of the Disks of the Simulated Spiral GalcLxies 

So far we have investigated the resulting morphology of our galaxy models via A^-body 
simulations. However, this procedure is somewhat incomplete because it is necessary to 
recover some observational counterparts from the simulated systems in order to test the 
reliability of our models. 

One of these counterparts is just the disk density profiles. In spite of numerical noise, 
very common in A^-body simulations, one expects that, in the cases where secular equilibrium 
plays an important role, the simulated systems maintain their density profiles, even when 
the mixing in the phase space occurs. 

Then, we calculate from the final snapshots (t=l Gyr), the density profiles of the disks 
in all cases presented in Sections 3.2-3.4. In Figure [TBI we display the final density profiles 
of the simulated disks and compare them with the initial profile. The method used here 
considers that the disk is composed by concentric rings in the xy-plane, where we count the 
A^* particles in each annulus and apply the formula p = N*/A, where A is the ring's area. 

Figure [TS] shows that our Newtonian model is stable, even for 1 Gyr of simulated time, 
since the initial and final density profiles are almost the same and present an exponential 
shape. 

On the other hand, spiral galaxies modeled with the YGP preserve the exponential disk 
profile only for A > 100 kpc. This lower limit for A w as also obtained in our previous work 



with elliptical systems f lBrandao fc de Araujd 1201 Obi ) 



For A < 100 kpc there is not an exponential density profile at the end of the simulation. 
For example, for A = 1 kpc the exponential d isk is completely destroyed. This result is 



consistent with the semi-analytical approach by Ide Araujo fc Mirandal (120071 ). 
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It is worth noting that the density profiles for the high- and low-resolution simulations 
are similar, that is why we show in Figure [18] only the results for the low resolution. 



From our simul ations, it is possible to conclud e that the graviton mass is nig <^ 10 
g, in agreement with iBrandao fc de Araujd (j2010bl ). 
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The density profile diagnosis then helps one to show how to constrain alternative 
gravitational potentials, since not necessarily the simulated images of a spiral galaxy and 
nor its rotation curve allow to constrain. 



4. CONCLUSIONS 

We give a recipe to probe alternative theories of gravitation in the non-relativistic 
regime, using A^-body simulations to model spiral galaxies. As an example, we use the 
recipe given here for the YGP. It is worth stressing that this recipe can be applied to any 
other alternative theory of gravitation in which the superposition principle is valid. In fact, 
in forthcoming studies we will probe other alternative potentials using the recipe given here. 

Basically, one has just to modify the code where the gravitational potential and the 
corresponding gravitational acceleration are taken into account. Then, one has to model 
a galaxy, which is initially made consistent with observations and run the simulation with 
the modified (non-Newtonian) code. For the alternative theory be reliable, the simulated 
galaxy at, say, 1 Gyr should resemble a true galaxy, i.e., it must have, for example, a disk 
consistent with that inferred from observations. Models presenting the destruction of galaxies 
or presenting bizarre morphologies would mean that the potential under consideration would 
be unreliable. 



Although Ide Araujo fc Mirandal (120071 ) studied, for example, disk galaxies under YGP, 
they used analytical arguments (centrifugal equilibrium), while our galaxies behave like 
"living" systems because they are composed by thousands of self-gravitating particles. In 
this way, this can be considered a reliable and strong test due to the fact that A^-B ody 



systems are very sensitive to chaos and complex phenomena (IBinney &: Tremainell2008l ). 



In this work, we have studied some models of late-type systems to probe the YGP and 
to constrain the Yukawian A parameter. As expected, we have seen that if A is much larger 
than the characteristic dimensions of spiral galaxies, the YGP and the Newtonian models 
yield the same results. On the other hand, for A smaller than the characteristic dimensions 
of galaxies, the YGP and the Newtonian models yield very different results. Moreover and 
more importantly, YGP galaxies for small values of A do not reproduce the rotations curves 
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and the surface density profiles observed in spiral galaxies. 

As a general conclusion, if YGP were reliable we should have A > 100 kpc, otherwise, we 
could not see late-type systems in the universe. This value of A is larger than that inferred 
from the solar system's constraints and it could be considered a good estimative, since with 
such a value the simulated galaxies remain "alive" for billions of years and look like their 
observational counterparts. 

Last but not the least, we intend in the near future to follow an interesting suggestion 
given by an anonymous referee, namely, instead of starting our simulations with a Newtonian 
model we will start with collapsing halos under the Yukawa potential to investigate the 
characteristics of the equilibrium halo profile, in particular if it resembles to some of the 
halo profiles discussed in the literature. 

C.S.S.B. and J. C.N. A. thank the Brazilian agencies CAPES and FAFESF for support. 
J. C.N. A. also thanks the Brazilian agency CNPq for partial support. C.S.S.B. also thanks 
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exchange of e-mail messages concerning the realization of disks. Finally, we thank the referee 
for the careful reading of the paper, the criticisms, and the very useful suggestions which 
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Fig. 1. — Rotation curves of our galaxy model in cylindrical coordinates. The abscissa 
shows the distance from the center, in kpc, in the plane of the disk. The ordinate shows the 
velocities for the disk, the halo, and the whole galaxy. 
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Fig. 2. — Top left: phase space for initial snapshot data. Top right: energy conservation of 
the simulation. Bottom left: phase space for final snapshot data at 1 Gyr. Bottom right: 
Rotation curves for initial (solid line) and final (dashed line) snapshots. 



- 19 - 




Fig. 3. — Newtonian disk at 2— projection at 0, 0.33, 0.66, and 1 Gyr of simulated time 
(indicated in the respective boxes). 
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Fig. 4. — First 320 Myr of simulated time to the Newtonian disk at 2;— projection. Time is 
indicated in the respective boxes. 
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Fig. 5. — Top left: phase space for the initial snapshot for the Yukawian disk simulation 
with A = 100 kpc. Top right: energy conservation of the simulation. Bottom left: phase 
space for final snapshot data at 1 Gyr. Bottom right: Rotation curves for the initial (solid 
line) and the final (dashed line) snapshots. 
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Fig. 6. — Disk at 2;— projection at 0, 0.33, 0.66, and 1 Gyr of simulated time (indicated in 
the respective boxes) for A = 100 kpc. 
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Fig. 7. — First 320 Myr of simulated time to the Yukawian disk at 2;— projection for A = 100 
kpc. Time is indicated in the respective boxes. 




Fig. 8. — Same as in Figure [5] for the high-resolution simulation for A = 100 kpc. 
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Fig. 9. — Same as in Figure [6] for the high-resolution simulation for A = 100 kpc. 
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Fig. 10. — Top left: phase space for the initial snapshot for the Yukawian disk simulation 
with A = 10 kpc. Top right: energy conservation of the simulation. Bottom left: phase 
space for final snapshot data at 1 Gyr. Also shown is a zoom of the final rotation curve, for 
comparison. Bottom right: Rotation curves for the initial (solid line) and the final (dashed 
line) snapshots. 
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Fig. 11. — Disk at 2— projection at 0, 0.33, 0.66, and 1 Gyr of simulated time (indicated in 
the respective boxes) for A = 10 kpc. 
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Fig. 12. — First 320 Myr of simulated time to the Yukawian disk at 2;— projection for A = 10 
kpc. Time is indicated in the respective boxes. 
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Fig. 13. — Same as in Figure [TO] for the high-resolution simulation for A = 10 kpc. 
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Fig. 14. — Same as in Figure [TT] for the high-resolution simulation for A = 10 kpc. 
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Fig. 15. — Top left: phase space for the initial snapshot for the Yukawian disk simulation 
with A = 1 kpc. Top right: energy conservation of the simulation. Bottom left: phase space 
for final snapshot data at 1 Gyr. Bottom right: rotation curves for the initial (solid line) 
and the final (dashed line) snapshots. Also shown is a zoom of the initial rotation curve for 
comparison. 
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Fig. 16. — Disk at 2— projection at 0, 0.33, 0.66, and 1 Gyr of simulated time (indicated in 
the respective boxes) for A = 1 kpc. 
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Fig. 17. — First 320 Myr of simulated time to the Yukawian disk at 2;— projection for A = 1 
kpc. Time is indicated in the respective boxes. 
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Fig. 18. — (I) Analytical exponential profile (solid line). Points represent the particle counts 
per unit area, obtained from the initial snapshot. The dotted line represents the radial profile 
of the final (at 1 Gyr) snapshot for the Newtonian run. (II) The initial profile (solid line), the 
final profile for A = 1 kpc run (dotted line), and the final profile for the A = 10 kpc (dashed 
line). (Ill) The initial profile (solid line) and the final profile for A = 100 kpc (dotted line). 
(IV) The initial profile (solid line) and the final profile for A = 1000 kpc (dotted line). 



